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A two-dimensional dipolar Fermi gas in harmonic trap under rotation is studied by solving ” ab ini¬ 
tio” Kohn-Sham equations. The physical parameters used match those of ultracold gas of fermionic 
molecules, a prototype system of strongly interacting dipolar quantum matter, which has 
been created very recently. We find that, as the critical rotational frequency is approached and the 
system collapses into the lowest Landau level, an array of tightly packed quantum vortices develops, 
in spite of the non-superfluid character of the system. In this state the system looses axial symmetry, 
and the fermionic cloud boundaries assume an almost perfect square shape. At higher values of the 
filling factor the vortex lattice disappears, while the system still exhibits square-shaped boundaries. 
At lower values of the filling factor the fermions become instead localized in a ’’Wigner cluster” 
structure. 

PACS numbers: 


Recent developments in the field of ultracold dipo¬ 
lar atoms have boosted investigations of many-body ef¬ 
fects associated with long-range interactions (see, e.g., 
the review in Ref. Q)- Recently, the creation of ultra¬ 
cold dipolar gas of fermionic molecules with large in¬ 
trinsic dipole moments has been achieved [3, 0 , open¬ 
ing the way to explore the intriguing many-body physics 
of correlated Fermi systems associated with the long- 
range, anisotropic nature of dipolar interaction between 
molecules 0, which include topological superfluidity d, 
[^, interlayer pairing between two dimensional systems 
and the formation of dipolar quantum crystals^]. 

Two-dimensional (2D) dipolar systems are particularly 
interesting, since the lifetime of heteronuclear Feshbach 
molecules with permanent electric dipole moment in¬ 
creases by the confinement in two dimensions Q. Indeed 
such polar molecules can have very large dipole moments, 
of the order of 1 Debye, allowing to access the regime of 
strong correlations in a controllable way. A 2D dipolar 
Fermi liquid, which is stable at low density, is expected 
to convert into a Wigner crystal at high densities P, [^. 
For intermediate values of the interaction strengths an 
instability at finite wave vector is predicted, driving the 
system to a ’’stripe”- phase [lol - [l^ (see also Ref.[lj and 
references therein). Recent Quantum Monte Carlo study 
confirmed the liquid-solid transition at high coupling but 
found that the stripe phase is never energetically favored 
0 . 

A particular fascinating route towards the realization 
of strongly correlated system of ultra-cold gases (either 
bosonic or fermionic) is the use of rapidly rotating har¬ 
monic traps. When the rotational frequency approaches 
the trap frequency, i.e. just below the limit of centrifugal 
instability, the single particle energy spectrum becomes 
highly degenerate and hence the kinetic energy of the 
Fermi system is much reduced, thus enhancing the role 
of the interparticle interactions. 


A uniform rotation with angular velocity approach¬ 
ing the centrifugal limit is in fact formally equivalent 
to a magnetic field (in the rotational frame) that re¬ 
group single-particle states into discrete, highly degen¬ 
erate Landau Levels (LL). Such equivalence, which holds 
in 2D and in the presence of an harmonic trapping poten¬ 
tial only, is embodied in the following formal identity [l3| 
involving the many-body Hamiltonian of the (interact¬ 
ing) system in the rotating reference frame: 

N 

= X nf + ^{ujI - n^)r^] + V (1) 

i=l ^ 

where V is the interaction energy, fl is the rotation fre¬ 
quency and Li 2 is the projection of the angular mo¬ 
mentum of the i-th particle along the z axis. Here 
Vi = XiGx + UiGy is the position vector of the z-th 
particle. When Fl = coh the non-interacting part re¬ 
duces to the Landau Hamiltonian of particles with mass 
M and charge e moving in a constant magnetic field 
B = Bbz, of strength B = 2MFl/e. The eigenvectors 
of the non-interacting part span Landau levels with en¬ 
ergies Cn = tiu!c{n + 1/2), where Wc = 2D. 

We consider in the following a two-dimensional, spin- 
polarized dipolar Fermi gas, characterized by an inter¬ 
action term in Eq.(l) V = Here d is the 

magnetic dipole moment of an atom/molecule and r, r' 
are coordinates in the 2D x — y plane. Being the dipole 
moments aligned parallel to the z-axis, the (long-range) 
pair potential is purely repulsive. The range of the 
dipole-dipole interaction is characterized by the length 
ro = Md?/h^. 

Exotic forms of vortex lattices, e.g. square, stripe- 
and bubble-” crystal” lattices are expected in rotating 
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Bose Einstein condensates when the critical rotational 
frequency is approached 15-171. Rotating dipolar Fermi 
gases have been proposed TMlii as suitable candidates 
to realize Laughlin-like state and more exotic quantum 
liquids, as well as their crossover behavior to Wigner 
crystals. At variance with the case of a non-rotating 
dipolar Fermi gas in a 2D trap, where the crystalline 
state becomes energetically favored at high densities, in 
the case of a fast rotating dipolar gas, the situation is 
reverse dfl: rotating dipoles in the Lowest LL (LLL) be¬ 
have similarly to electrons, where the crystalline phase is 
stable at low densities. Indeed, it has been shown [ 2 ^ 1 ^ 
that a rapidly rotating polarized 2D dipolar Fermi gas 
undergoes a transition to a crystalline state, similar to 
the two-dimensional Wigner electron crystal in a mag¬ 
netic field, for sufficiently low value of the filling factor, 
V < 1/7. Here v = 2TTpnp (where np is the areal density 
of the fermionic system and I = y^H/Muih is the magnetic 
length) gives the fraction of the occupied LLL. At filling 
factor p = 1/3 the system is instead well described in 
terms of fractional quantum Hall-like states[l6j|. 

Density Functional Theory (DFT), which is perhaps 
the most widely used and successful technique in elec¬ 
tronic structure calculations of condensed matter sys¬ 
tems, has only recently entered the field of cold gases 
as a useful computational tool which goes beyond the 
mean-field description by taking into account correlation 
effects, and thus is capable to yield quite accurate results 
in agreement with more microscopic (but also much more 
computationally expensive) ^proaches. The well-known 
Kohn-Sham (KS) mapping[22[ of the many-body prob¬ 
lem into a non-interacting one make this approach appli¬ 
cable in practice, often within the so-called Local Density 
Approximation (LDA)[2^. Recently, KS-DFT has been 
applied to cold atomic Fermi gases in optical lattices [ 2 ^ 
and to the study of unitary trapped Bose gas [l^. DFT 
approaches have been used recently to describe a Fermi 
dipolar system in various ’’single-orbital” a ppr oximations 
(Thomas-Fermi [^, Thomas-Fermi-Dirac 251, Thomas- 
Fermi-von Weizsacker 2^, 271). In Ref. [ 2 ^ a parameter- 
dependent DFT-LDA approach was used to study small 
number of harmonically trapped fermions. A somewhat 
different density functional formalism, whose applicabil¬ 
ity is however limited to a small number of particles, and 
which is based on the self-consistent combination of the 
weak and the strong coupling limits, has been proposed 
to study the ground-state properties of strongly corre¬ 
lated dipolar and ionic ultracold bosonic and fermionic 
gases [^. 

Here we use the conventional KS approach, based on 
accurate description for the correlation energy of the 
dipolar system as provided by Diffusion Monte Carlo 
calculations [l^. Our approach does not require any ad¬ 
justable parameter, and thus belongs to the family of 
the ”ab initio” methods well known in the electronic 
structure community. The Kohn-Sham formulation 


of Density Functional Theory for an inhomogeneous 
system of N interacting particles with mass M is based 
on the following energy functional of the density which 
includes the exact kinetic energy of a fictitious non¬ 
interacting system and the interaction energy functional 
Ehfc- 


EksIp] = -^ E / </>:( r ) V ^(/.( r)dr + Ehpc[p] ( 2 ) 

The i = 1, N} are single-particle orbitals, form¬ 

ing an orthonormal set, (/'il/'j) = Sij, filled up to the 
Fermi level. The total density of the system is p(r) = 

Ehfc is the sum of the direct-l-exchange dipolar in¬ 
teraction term (usually termed ” Hartree-Fock” energy, 
Ehf) and the correlation energy {Ec)- The Hartree- 
Fock energy of a dipolar Fermi gas in 2D has two contri¬ 
bution. The first, for the homogeneous system of surface 
density p = N/A, is: 




( 3 ) 


The second term is non-local in nature and is given 
by[2i[26| 


4'f = I drp{r) J dv' j ^fce-*'^ (---')p(r') 

( 4 ) 

This term vanishes in the uniform limit, while the neg¬ 
ative sign crucially lower the total energy of the system 
in inhomogeneous configurations. This term has been 
shown to be essential to stabilize structures such as one¬ 
dimensional stripe phases and the Wigner crystal that is 
expected at high densities 

In the following we will treat and Ec within the 
Local Density Approximation (LDA), i.e. 


E% + Ec = j [^d^y/^p{rf/'^ -b p{r)ec{p{r))]dr (5) 

where cc{p) is the correlation energy per particle of the 
homogeneous system of density p, as obtained from the 
(virtually exact) Diffusion Monte Carlo calculations of 

Ref. [13 . 

The total energy functional in the co-rotating frame 
with constant angular velocity H (where the dipolar sys¬ 
tem appears at rest) and in the presence of an isotropic 
harmonic trapping potential of frequency ujh, U{r) = 
+ y^)) is given by: 


E[p] = Eks[p] + J drp{r)U(r) - n{L^) ( 6 ) 


















Here (L^) is the total angular momentum of the system. 
Constrained minimization of the above functional leads 
to the coupled KS eigenvalues equations 
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2M 


VKs](t>i(,i^) = e*(/'n(r) 


(7) 


where 

VKs{r) = ec(p(r)) + 

op 9 

-2nd^ Jdr'J ( 8 ) 

and Lz = —ih{xd/dy — yd/dx). 

We seek for stationary solutions {4>i{r),i = l,iV} 
by pro pag ating in imaginary time the time-dependent 
version(2^ of the KS equations 0 . Both the density 
and the orbitals (j)i have been discretized in cartesian co¬ 
ordinates using a spatial grid fine enough to guarantee 
well converged values of the total energy. The orthogo¬ 
nality between different orbitals has been enforced by a 
Gram-Schmidt process. The spatial derivatives entering 
Eq.© have been calculated with accurate 13-point for¬ 
mulas, while Fast-Fourier techniques have been used to 
efficiently calculate the non-local term entering the KS 
potential Vks- 

We take in our calculations d = 0.8 Debye, which is 
appropriate to 023 molecules in the experimental 

realization of Ref. 3 . The mass is that of a K 4 oNa 23 
molecule. The range of the potential is thus ro = 
Md?/h^ ^ 0.6 ym ^ 0.2 an, an = \/h/2Mwh being the 
oscillator length. In the ground-state of the non-rotating 
system, the adimensional interaction strength character¬ 
izing the system is kpr^ ^ 0.9 (where kp = \/'knpmax is 
the Fermi wavevector of the 2D system at a density equal 
to the maximum density in the center of the trap), i.e. a 
relatively weak value which can easily be achieved in ex¬ 
periments. The interparticle distance (r) is larger than 
the range of the interaction, being (r)/ro 3.6. The 
corresponding dipolar interaction energy Ed = d'^/(r)^ 
approaches 20% of the local Fermi energy h^kp/2M. In 
spite of the relatively weak coupling, as a consequence 
of the rotation, strong correlation effects will show up 
in the density distribution of the calculated stationary 
states, as shown in the following. We consider systems 
with up to fV = 200 fermions. 

We show in Fig U the evolution of the calculated single 
particle KS eigenvalues Ci for the case N = 100, as the 
rotational frequency H approaches from below the har¬ 
monic frequency ujh. At H = 0.999 it appears that all 
the energy levels collapse into a single level, the (highly 
degenerate) LLL. 

It is instructive to follow how the density of the system 
evolves as H is increased. This is shown in Fig© where 
the densities of selected configurations corresponding to 
different values of D are displayed. 



FIG. 1: Calculated KS eigenvalues for N = 100 fermions, for 
different values of the rotational frequency D 



FIG. 2: (Colors online) Density of the dipolar system {N — 
100) at selected values of Q,/ujh. From left to right and from 
top to bottom: Q/tOh = 0.96, 0.98, 0.99, 0.999. The x-axis 

coordinates are in units of an, while the density is in units of 
-2 

“h ■ 


For low values of D (D = 0 included) the calculated 
stationary states have circular symmetry, and the density 
has the familiar, almost featureless shape of a trapped 
cold gas cloud. As soon as D approaches uJh, however, 
a ring of equally spaced deep dimples develop close to 
the periphery of the cloud, while the systems looses its 
axi-symmetric shape. Eventually, very close to lUh, an 
array of tightly packed vortices develops, similarly to the 
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Abrikosov lattice of vortices in rotating superfluids, while 
the system boundaries acquire a surprising square shape. 

The calculated current density in the state with 
fl/uJh = 0.999 in Figl2] appears indeed to be circulat¬ 
ing around the zero-density minima (black dots in the 
last panel of FigI5]), as expected for a vortex array. The 
total angular momentum {L^) shows also the typical be¬ 
havior associated to the nucleation of quantum vortices, 
i.e. a sequence of rounded steps (with amplitudes ^ Nh) 
with increasing rotation frequency as more vortices 
are nucleated in the system during the minimization pro¬ 
cess leading to the stationary state shown in the last 
panel of Figl^ The average distance dy between vor¬ 
tices in the structure shown in Figl2]) is ^ 20% larger 
than the one calculated (assuming a triangular vortex 
lattice of areal density riy) using Feynman’s formula 13ll. 
ny = 2/V^dl=]Vm/Trh. 


Quantized vortices in Fermionic cold gases are usually 
associated to pairing interactions, as in the BCS side of a 
unitary Fermi gas(32|, where they are considered the hall¬ 
mark of the superfluid character of the system. The pres¬ 
ence of vortices in a system with purely repulsive interac¬ 
tions, like the one studied here, has been predicted to oc¬ 
cur in fermion systems with purely repulsive interaction 
such as quantum dots, where the rotation is induced by 
an external magnetic field (see for instance RefIsSl. In¬ 
direct evidence of vortices in ultra small fermion droplets 
{N = 6) with aligned dipoles has been provided in Refl34 . 
However, due to the implicit symmetry constraints in the 
calculations of Refjs^ multi-vortex structures like the one 
shown in Fig l2] did not show up in the calculated density 
profiles. 

A striking feature of the N = 100 system in the LLL 
is the lack of axial symmetry represented by the un¬ 
usual square-shaped boundaries. This seems to be inti¬ 
mately connected with the interactions between fermions: 
the system shown in Fig 121 under the same conditions 
but with no interactions between fermions, exhibits den¬ 
sity profiles with circular symmetry all the way up to 
ujh- Deviations from axi-symmetric configurations in 
isotropic trapping have been found in fast rotating BEG 
at overcritical rotation (^. as a consequence of the in¬ 
teratomic forces. Stable, non axi-symmetric multi-lobed 
shapes also characterize the fast rotation of classical liq¬ 
uid droplets [1^. 

The configuration shown in the lowest panel of Figl2] 
corresponds to a filling factor v ~ 0.77. By decreasing 
the number of fermions in the trap we can reach lower val¬ 
ues of zz. One example is shown in Figl3l where N = 13. 
Again, as the centrifugal limit is approached, station¬ 
ary configuration with an increasing number of vortices 
are found: vortices enter the fermion droplet from the 
low density periphery (a mechanism common to BEClai 
and Helium-4 [38[). AsVl ^ coh (last panel in Figl3]), how¬ 
ever, a completely different pattern shows up, resembling 
a cluster of localized particles (albeit with a partially 



FIG. 3: (Colors online) Density of the dipolar system {N = 
13) at selected values of Q/uJh From left to right and from 
top to bottom: id/oJk = 0.90, 0.96, 0.98, 0.999. The x-axis 
coordinates are in units of an, while the density is in units of 
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FIG. 4: (Golors online) Density of the dipolar system at 
Q/ojh = 0.999 for N = 160 (left panel) and N = 200 (right 
panel), corresponding to filling factors zz = 1.08 and zz = 1.25, 
respectively. 


melted second shell). This configuration is characterized 
by zz ^ 0.18. We take this as a clear evidence of forma¬ 
tion of a ’’Wigner cluster” structure for sufficiently low 
values of the hlling factor. 

Higher values of zz can be conversely achieved by in¬ 
creasing the fermions number. In this case the vortex 
lattice disappears, and the smoother structures shown in 
Fig 0] develop. Here, zz = 1.08 and zz = 1.25, respectively 
(corresponding to IV = 160 and N = 200 fermions). Note 
however that the peculiar square-shaped boundaries re¬ 
main even at higher values of the filling factor. 

Although small numbers of cold trapped atoms, like 
the ones considered here, can nowadays be achieved 
in experiments [ 3 ^, we expect that the surprising phe¬ 
nomenology revealed by our calculations should be 
present also in larger systems of rotating dipolar 
fermionic molecules in quasi 2D harmonic traps. Achiev- 
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ing rotation frequencies ft ^ 0.999 w/j is a challeng¬ 
ing task, but definitely within the reach of current 
experiments [l^. Due to the relatively high contrast of 
the vortex array shown in Figl^l its observation should 
be possible by direct imaging of the atomic cloud after 
expansion (the vortex lattice in quantum dots still awaits 
experimental detection). 

The author thanks S.Giorgini and N.Matveeva for hav¬ 
ing shared their DMC numerical results, and S.Giorgini, 
L.Salasnich, F.Toigo, J.Boronat and M.W.Cole for useful 
discussions and comments. 
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